test_mod1a <- lm(java_index_skd_bw_1 ~ fail_skd, data = data_1 %>% filter(abs(forcing_skd) < 5))
test_mod1b <- lm(java_index_skd_bw_1 ~ fail_skd + AGE + factor(GENDER) + factor(urb_rur) + factor(islam), data = data_1 %>% filter(abs(forcing_skd) < 5))

test_mod2a <- lm(non_java_index_skd_bw_1 ~ fail_skd, data = data_1 %>% filter(abs(forcing_skd) < 5))
test_mod2b <- lm(non_java_index_skd_bw_1 ~ fail_skd + AGE + factor(GENDER) + factor(urb_rur) + factor(islam), data = data_1 %>% filter(abs(forcing_skd) < 5))

test_mod3a <- lm(regional_index_skd_bw_1 ~ fail_skd, data = data_1 %>% filter(abs(forcing_skd) < 5))
test_mod3b <- lm(regional_index_skd_bw_1 ~ fail_skd + AGE + factor(GENDER) + factor(urb_rur) + factor(islam), data = data_1 %>% filter(abs(forcing_skd) < 5))

test_mod4a <- lm(religious_index_skd_bw_1 ~ fail_skd, data = data_1 %>% filter(abs(forcing_skd) < 5))
test_mod4b <- lm(religious_index_skd_bw_1 ~ fail_skd + AGE + factor(GENDER) + factor(urb_rur) + factor(islam), data = data_1 %>% filter(abs(forcing_skd) < 5))

test_mod5a <- lm(corruption_index_skd_bw_1 ~ fail_skd, data = data_1 %>% filter(abs(forcing_skd) < 5))
test_mod5b <- lm(corruption_index_skd_bw_1 ~ fail_skd + AGE + factor(GENDER) + factor(urb_rur) + factor(islam), data = data_1 %>% filter(abs(forcing_skd) < 5))

test_mod6a <- lm(national_index_skd_bw_1 ~ fail_skd, data = data_1 %>% filter(abs(forcing_skd) < 5))
test_mod6b <- lm(national_index_skd_bw_1 ~ fail_skd +  AGE + factor(GENDER) + factor(urb_rur) + factor(islam), data = data_1 %>% filter(abs(forcing_skd) < 5))

observations <- c(nobs(test_mod1a), nobs(test_mod1b), nobs(test_mod2a), nobs(test_mod2b), 
                  nobs(test_mod3a), nobs(test_mod3b), nobs(test_mod4a), nobs(test_mod4b),
                  nobs(test_mod5a), nobs(test_mod5b), nobs(test_mod6a), nobs(test_mod6b))

test_mod1a <- coeftest(test_mod1a, vcov=vcovHC(test_mod1a,type="HC0"))
test_mod1b <- coeftest(test_mod1b, vcov=vcovHC(test_mod1b,type="HC0"))

test_mod2a <- coeftest(test_mod2a, vcov=vcovHC(test_mod2a,type="HC0"))
test_mod2b <- coeftest(test_mod2b, vcov=vcovHC(test_mod2b,type="HC0"))

test_mod3a <- coeftest(test_mod3a, vcov=vcovHC(test_mod3a,type="HC0"))
test_mod3b <- coeftest(test_mod3b, vcov=vcovHC(test_mod3b,type="HC0"))

test_mod4a <- coeftest(test_mod4a, vcov=vcovHC(test_mod4a,type="HC0"))
test_mod4b <- coeftest(test_mod4b, vcov=vcovHC(test_mod4b,type="HC0"))

test_mod5a <- coeftest(test_mod5a, vcov=vcovHC(test_mod5a,type="HC0"))
test_mod5b <- coeftest(test_mod5b, vcov=vcovHC(test_mod5b,type="HC0"))

test_mod6a <- coeftest(test_mod6a, vcov=vcovHC(test_mod6a,type="HC0"))
test_mod6b <- coeftest(test_mod6b, vcov=vcovHC(test_mod6b,type="HC0"))


table <- list(test_mod1a, test_mod1b, test_mod2a, test_mod2b, test_mod3a, test_mod3b, 
              test_mod4a, test_mod4b, test_mod5a, test_mod5b, test_mod6a, test_mod6b)

note_text <- paste(" Beta coefficients from OLS regression with and without controls. Standard errors were calculated using the Huber-White (HC0) correction. 
                   The outcomes measure are indexed values capturing (1) Javanese preferentialism among Javans and (2) among non-Javans, (3) regional preferentialism,
                   (4) religious resentment, (5) perceptions of corruption, (6) national identification.")

table = stargazer(table, 
                  type = 'latex', 
                  title = "The Effect of Basic Competence Examination (SKD) Failure, Controls",
                  label = 'tab:testing_effect_controls',
                  model.names = F,
                  model.numbers = T,
                  digits = 3,
                  column.separate = c(4, 2, 2, 2, 2),
                  column.labels = c("Java. Pref.", "Reg. Pref.", "Relg. Resent.", "Corrup. Percep.", "Natl. ID"),
                  
                  multicolumn = T,
                  dep.var.labels = NULL, 
                  add.lines = list(c("Subset", "Javan", "Javan", "non-Javan", "non-Javan", "---", "---", "---", "---", "---", "---", "---", "---"),
                                   c("Observations", observations),
                                   c('Bandwidth', rep(c('1\\%'), 12))),
                  covariate.labels = c("Failed SKD", "Age", "Woman", "Urban", "Muslim"),
                  #star.cutoffs = c(0.05, 0.01),
                  #float.env = 'sidewaystable',
                  keep.stat = c("n"),
                  notes = NULL,
                  notes.align = 'l')


write_latex(table[-c(10, 11, 12, 18, 21, 24, 27, 30, 33, 40)], note_text, './_4_outputs/tables/table_a9.tex', .8)

